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We derive an analytical expression for extracting the gravitational waveforms at null infinity using 
the Weyl scalar ^4 measured at a finite radius. Our expression is based on a series solution in orders 
of 1/r to the equations for gravitational perturbations about a spinning black hole. We compute 
this expression to order 1/r 2 and include the spin parameter a of the Kerr background. We test the 
accuracy of this extraction procedure by measuring the waveform for a merging black-hole binary 
at ten different extraction radii (in the range r/M = 75 — 190) and for three different resolutions 
in the convergence regime. We find that the extraction formula provides a set of values for the 
radiated energy and momenta that at finite extraction radii converges towards the expected values 
with increasing resolution, which is not the case for the ‘raw’ waveform at finite radius. We also 
examine the phase and amplitude errors in the waveform as a function of observer location and 
again observe the benefits of using our extraction formula. The leading corrections to the phase 
are (7(1/r) and to the amplitude are (7(1/r 2 ). This method provides a simple and practical way of 
estimating the waveform at infinity, and may be especially useful for scenarios such as well separated 
binaries, where the radiation zone is far from the sources, that would otherwise require extended 
simulation grids in order to extrapolate the ‘raw’ waveform to infinity. Thus this method saves 
important computational resources and provides an estimate of errors. 
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I. INTRODUCTION 


Perturbation theory about black-hole backgrounds and fully nonlinear numerical simulations of the Einstein field 
equations provide complementary approaches to solving important problems in Relativity. A few examples of the 
synergy created by using the two together include the use of perturbative boundary techniques for a fully nonlinear 
simulation far from the sources as a way of propagating most of the radiation out the simulation domain [T] and the 
now classical Lazarus approach (2H3 that extracted spatial information from a short lived full numerical evolution to 
provide initial data for a subsequent perturbative evolution of a single rotating black hole. 

With the breakthroughs in numerical relativity EHE], complete simulations of inspiraling black-hole binaries became 
possible. However, even in this case, the spacetime far from the sources (more precisely, in the wavezone) can be 
described by black hole perturbation theory. Here we will exploit this fact to analytically propagate the waveform 
from a fully nonlinear simulation, but measured at finite distance from the sources, to null infinity. 


Over the last few years several of these waveform extraction techniques have been developed. The most straightfor¬ 
ward strategy would be to have the numerical domain extend very far from the sources and extrapolate the waveform 
measured at far distances to infinity. This can be achieved at reasonable computational efficiency using pseudo-spectral 
decomposition of the fields Pi, or by using multi-patch techniques mm- 

A more sophisticated waveform extraction technique, and one that produces a true gauge invariant signal, is Cauchy- 
Characteristic extraction (CCE) jT5ffT7] . In this technique, the metric and its derivatives on a timelike worldtube are 
used as inner boundary data for a subsequent characteristic evolution. As the characteristic evolution includes null 
infinity, the waveform obtained is exact (up to truncation error). A complementary approach to CCE is to evolve 
the spacetime on surfaces that are spacelike in the interior but asymptote to null slices that intersect, null 

infinity Bang. 


An alternative extrapolation method consists of using the results of perturbation theory to propagate waveforms 
obtained at finite radii (but in the radiation zone) to infinity. Treating the background spacetime as a perturbation of 
Schwarzschild, which will be accurate in the wavezone, leads to a simple explicit formula relating 7^4 at infinity with 
the finite radius ^4 and its time integral. For more details, see Ref. [20], Eq. (53). This method has been proven to 
correct for the next-to-leading 1/R 0 bs term in R 0 Robs) OH] using only a single observer radius and displays 
a significantly reduced level of extrapolation noise, when compared to the standard polynomial extrapolation. The 
errors produced by this method can be estimated by applying it to different extraction radii. We applied this method to 
the q = 10 case in Ref. [22] and found good agreement (but with significantly reduced noise) between the perturbative 
and standard extrapolation technique used in this paper. 


In this paper we expand upon this method by including higher-order [0(l/r 2 )\ and rotation effects for extracting or 
extrapolating fields from an intermediate distance to infinity via a perturbative, analytic expansion. This method is 
relatively simple to implement, yet it is accurate enough for most applications. This includes cases of large-separation 
binaries with long orbital periods leading to wave zones extending beyond several thousand M (see for instance 
Ref. [23] where evolution of a binary separated by 100 M led to waveforms with 6400M periods), as well as studies of 
the scattering of two black holes starting far apart to measure small scattering angles [24] and high energy collision 
of black holes [25], which require extractions at large distances from the sources. Another circumstance when this 
extraction method can be of use is when more physical scales need to be resolved. Such is the case when matter 
surrounds black holes [26] or when one tries to simulate a hybrid systems involving neutron stars and black holes or 
binary neutron stars m 


The paper is organized as follows. Section [TT] discusses the extraction of gravitational waves propaga ting as a 


perturbation on the asymptoti c Sch warzschild background with 0(1/r) corrections included in Sub-Sec. IIA and 


0(1/r 2 ) corrections in Sub-Sec. IIB In S ec.|lll| we include the effects of the spin in the background. Linear corrections 
in the spin in Sub-Sec. Ill A In Sub-Sec. IIIB we correct the extractions of the Weyl scalar ^4 for a nonconventional 
choice of the tetrad used in full numerical simulations. While in Sub-Sec. |III C| we collect together a definitive formula 
to include all effects together. This formula is capable of extracting numerical ^4 relatively close to the sources and 
extrapolate waveforms to infinity with accuracy, particularly for its phase and amplitude. Section [TV| contains explicit 
expressions for the radiated energy and momenta (along the z-axis) based on the extrapolated waveforms. In Sec. [V] 
we apply those equations into a case study of full numerical evolution of binary black holes. We choose ten extraction 
radii in the intermediate radiation region and evolve with three different resolutions in the convergence regime to 
study the effects of finite resolution on the extrapolated quantities. We finish the paper with a brief discussion in 
Sec. |VI| of the range of applicability of our results. 
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II. PERTURBATION IN A NONSPINNING BACKGROUND 

In this section we derive expressions relating the Regge-Wheeler-Zerilli (RWZ) [28, 29] functions at finite radius 
to their values on in the Schwarzschild (mass M ) black hole perturbation. Then, using these expressions, we 
derive expressions for the Weyl scalar ^4 at based on its values at finite radii. We always work in the first order 
perturbative regime, i.e. no quadratic terms in the perturbations around the black hole background are included, and 
expand the solutions asymptotically in powers of 1 /r. 


A. First-order corrections: (l/r)-terms 


The Weyl scalar, ^ 4 , in an asymptotically flat tetrad, like Kinnersley’s |50j, is related to the strain at large radii by 

lim r 04 = lim r(/i + — ih x ). (1) 


r—>• 00 r—> 00 

Similarly, the RWZ even and odd parity functions are related to the strain on at large radii by 


im 


— 2 Wm 5 


( 2 ) 


where and are the even and odd parity wave functions, respectively, and _ 2 bm denotes the spin(— 2 )- 

weighted spherical harmonics (see, e.g., review papers [3TH34] ). 

The asymptotic values of ^4 and the RWZ wavefunctions can be related to their values at finite radii by examining 
the asymptotic behavior of the relevant wave equations. For the RWZ wave equations we get, 


^ven/odd) = ^(even/odd)^ _ ^ + J (f - r *) + Q(l/r 2 ) , 


(3) 


for general £ modes, where H is the strain observed at infinity, and r* = r + 2 Mln[r/( 2 M) — 1 ]. An error due to 
finite extraction radii arises from the integral term in Eq. <i- Inverting the above relation, we have m 




(even/odd) 

im 


_ (even/odd) (± £(£ ~\~ 1 ) 


im 


/ ^L Ven/ ° dd) (^d + 0 (l/r 2 ). 


Similarly, if the Weyl scalar, 


f4 = ^<-2^, 


(4) 


(5) 


im 


satisfies the Teukolsky equation [35] in the Schwarzschild background spacetime, then the asymptotic behavior of 
is given by 


rVr = - r*) + - r*) + 0(l/r 2 ), 


(6) 


where the difference between H£ m and H£ m = H^ en ^ — iH^ d ^ defined from Eq. ^ is only a numerical factor, and 
we have the relation by using Eqs. © and © as 


* _ 

-H-im — 0 -H-irn • 


Inverting Eq. (| 6 |, we get 


rV4 m L = 


— 1)(^ + 2) 

2 r 


j dt[r^{ m (t,r)]+0(l/r 2 ). 


To see the phase and amplitude corrections by using the above formula, we assume 

H (even/°dd) ^ = ^ e xp(-M^ m (t _ r *)) , 

in Eq. ©• Then, the RWZ functions at a finite extraction radius are given by [36 1 


’bm en/0dd) = A im 


i + i (W + V 

2u lm r 


0(l/r 4 ) 


exp(-iui im (t ~ r*)) exp (i6</>e m ) + 0{l/r 2 ), 


(7) 

(8) 

(9) 

( 10 ) 
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where is defined as 

sin §4> im = ^ ^ + 0(1/r 2 ). (11) 

2 u>e m r 

Therefore, the phase correction from the perturbative formula has 0(l/r). On the other hand, from Eq. ( p~0| ) the 
amplitude correction will be 0(1/r 2 ) which we have ignored here. This result is consistent with Refs. [371l38| . and 
also has been observed in the black hole perturbation approach mm- This above analysis is also applicable to the 
Weyl scalar. In the next subsection, we extend the perturbative formula to order 1/r 2 . 


B. Second order corrections: (l/r 2 )-terms 

In this subsection, we discuss the next order correction of vpA™ in the 1/r-expansion, first on a Schwarzschild 


background. The starting point is the RWZ formalism and Eq. (|3j) is extended to order 1/r 2 . For the even parity 
function, we have 


dtH. 


(even) 


(*-0 + 


* ( z en) = ^z n \t-n + e -^ I 

3(£ 2 + £ + 2) M j dtH (even) ^ + 0(1/r 3 } , 


{£-!)£{£+!){£ +2) 


8 r 2 


/ JdtdtHj,Z n \t~r*) 


2 (£ - 1)(£ + 2) r 2 
and for the odd parity function, 


( 12 ) 


*L dd) = - r*) + J dtH^ d \t - r*) + 


(l-l)l(l+l)(l + 2) 
8r 2 


3 M 
2 r 2 


/ 


J JdtdtHj£ d \t-r*) 


dtH, 


(odd) 

£m 


(t - r*) + 0(l/r 3 ). 


(13) 


There is a difference between the even and odd parity functions at order 1/r 2 due to the difference in the potentials 
of the RWZ equations. 

Next, we convert the above even and odd parity functions into the Weyl scalar. Using Eqs. (C.l) and (C.2) in 
Ref. [33] and taking care of the definitions in Eqs. 0 and ([2]), we obtain 


V> 4 to = VII 


'{£ + 2 )! 


Um 2r \l (£ — 2)\ 


agrH* - o + - ™ +t) 6SrH' - r*> 


+ 


(£-!)£(£+ !)(£ +2) H(even) 


Sr 2 


£m 


3 M • (even) , 


\t - r*) + - r*) + 0( 1/r 3 ) 


^4 Zm = VI 


—i (£ + 2)1 


2 r I/ (£ — 2)! 


H, 

+ 


+ 1 // — ) gg J) («-r) 


(odd) 

£m V L 1 J 2 ±± £m 

(£_ — i)l(l + i)(l + 2) jj (odd) 


8r 2 


(^-^) + ^t dd) (^-^) + 0(iA 3 ) 

/,+/- 


(14) 


where the dot denotes the derivative with respect to the retarded time (t — r*). The functions Vb/m are defined in 
Eq. (13) of Ref. [33] and are the symmetric and antisymmetric Weyl scalar fields, respectively. It is natural to have 
the same asymptotic behavior for the Weyl scalar fields. 

Combining the above ^Xim as V4 m = + ^Um j we obtain the extension of Eq. |d) as 


ri ,em_fr (f „.(t-l)(£ + 2 )~ (^-l)^ + l)(* + 2) - 

^ U 4 — ^ j t 0 ^ ) \ Q H.£ rn [t ^ ) 


2 r 

+ y^H em (t - r*) + 0(l/r 3 ). 


8r 2 


(15) 


Inverting this equation, the perturbative formula extended to order 1/r 2 becomes 


rV4 m | 


, =r^ m (f,r)- 


(l-l)(l + 2) 
2 r 


(£-l)(£ + 2)(^ 2 +^-4) 

S^ 2 


J dt[riP{ m (t,r)} 

J J dtdt[rZ 4 m (t,r)}-^ J dt[ry 4 m (t,r)}+0(l/r 3 ). (16) 


The above relation is valid for the extrapolation of the ^4 i n the Kinnersley tetrad. Next we will consider the 
corrections due to spin and the use of a tetrad used in numerical relativity (NR) at a finite r and its decomposition 
into (£, m) modes. 
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III. IN SPINNING BACKGROUND 


In this section, we include the spin dependence in the Teukolsky formalism [35] of the Kerr (mass M and Kerr 
parameter a) black hole perturbation. It is noted that the wave function in the Teukolsky equation is _ 2 ^ = 
(r — iacos#) 4 ^. Here, we ignore 0(l/(u;r) 3 ) and 0((au:) 2 ) to derive a perturbative extrapolation formula from 
the frequency domain analysis. For example, when the extraction radius is r = 100M for a = M, the rough error 
estimation gives 1 /(u;r) 3 = 0.237% and (acj) 2 = 0.563% at uj = 0.075/M, respectively. This frequency is a reference 
to produce a hybrid post-Newtonian (PN)-NR waveform for the (£, m) = (2, 2) mode in the Numerical INJection 
Analysis (NINJA) project [41]. 


A. Background spin correction 


First, we focus on the Teukolsky’s wave function. 


2^ = / dtv^2- 

J tm 


2*w(r)- 2 «M)e 


—iuit 


(17) 


It is noted that we have used the spin-weighted spheroidal harmonics (_ 2 *S^(#,0)) in the Teukolsky formalism, 
while the spin-weighted spherical harmonics _ 2 are used in the NR simulations. The spin-weighted spheroidal 
harmonics, which are the solution of the angular Teukolsky equation, can be expanded as [42] 


- 2 S% = - 2 Yem + aujJ2 4m~*Y Vm + 0((au ) 2 ), 

where the coefficient c\ m has a non-zero value only for £' = £ ± 1, 


(18) 


= - — A 

u tm £2 ' 


'(£ + 2)(£ - 2)(£ + m){£ - m) £+1 _ 

’ C £m ~ 


' {£ + 3)(£ — 1)(£ + m + 1){£ — m + 1) 


(2£ — 1)(2£ + 1) ’ (^ + 1) 2 V (2£ + 1)(2£ + 3) 

The radial Teukolsky equation in the frequency domain gives the asymptotic solution, 


(19) 


— 2 


1 


+ iM + 


-4 ima i(£- 1)(^ + 2)\1 
J{£+1) + 2cu 
(£ + 2 ) (£- 1 ) 


cu£(£ + 1) 


ma + 


3 iM l£(£-l)(£ + 2)(£+l)\ 1 


2 to 8 


uj* 


3 +0(1/^, (own 


H, 


irrvjj 


i + 


-4 ima *(^-l)(^ + 2)\ 1 _ 1 1 {£ - 1) {£ + 2) {£ + 1) 1 
£ (£ + 1) 2 to Jr 8 to 2 r 2 


+ [higher order] 


H, 


tmuj • 


( 20 ) 


Hirnuj is related to the waveform at infinity. In the last line of the above equation, we have ignored various cross terms 
which are included in [higher order], (Muj)(cluj) /(cjr) 2 , ( auj)/(ujr ) 2 and (Mcj)/(u;r) 2 where we assumed that M and a 
are the same order. 

Inserting _ 2 ^mwW and _ 2 into Eq. (p~7|), we have 


2^ = 


/ 

tm ~ 


i + 


4 ima i(<-l)(f + 2)\ 1 1 £ [£ - 1) [£ + 2) [£ + 1) 1 


*(*+1) 


+ 


2w 




r 8 

—iut 


H-tmu-'lXt n 


[higher order] . 


( 21 ) 


Therefore, the spin-weighted spherical harmonic expansion becomes 
-2 


/ 


= / dui 


1 + 


4«mo i(^ —1)^ + 2)^ 1 U(^-l)(^ +2)^+1) 1 


^ (£ + 1) ' 2w 

+aw (c^ +lm ^+i mw -(- 


r 8 

—iu)t 


Hi 


imuj 


+ [higher order] . 


( 22 ) 
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B. Use of the full numerical tetrad 


Eq. (|8| in the nonspinning case relates the Weyl scalar -0 4 at a finite radii with the scalar at infinity. The preferred 
tetrad in perturbation theory of black holes is the Kinnersley tetrad [30] that make use of the algebraic specialty 
of the Kerr background where 7 /q vanishes. On the other hand, in full numerical relativity, the lack of a reference 
background makes this choice ambiguous and another tetrad, labeled ‘NR’, is conventionally used. This variant of 
the ‘psikadelia’ tetrad is described in Ref. [41. 

Using Eq. (2.15) in Ref. [7], we check the tetrad dependence. Assuming the peeling theorem (^4 = [r^J/r, ip 3 = 
[r 2 ^)/r 2 , ip 2 = [r 3 / 02 ]/ , r 3 , 'ipi = [ r4 Vh]/ r4 > V’o = [r^o]/?" 5 , where the functions in the square bracket are order r° for 
large r), we have 


r^f n 


= l K! - - lr *' 


ri 


1 a (7 a[r^ R ] cos 2 9 — 3 a[r^ R ]) 
4 


+i 


5 0 [rip. 




1 a (4 sin#[r 2 ?/> 3 fR 
4 


+ 1 


[rip, 


NRl 


COS 


9M)' 


£>(l/r 3 ). 


(23) 


After recasting the relationship between _ 2 ^ and ^4 in terms of the NR ^ 4 , we get 

_ 2 4> = ^ r 3 [r^ R ] — (M + iacosO) [r^ R ] r 2 + (2 iMa cos 0[np^ R ] — ia sin 0[r 2 ip^ R ]) r + 0(l/(o;r) 0 , (acu) 2 ) 

= ^ r 3 [r0^ R ] — (M + ia cos 0) [r^ R ] r 2 + [higher order] , (24) 

where we have ignored various cross terms again. The spin-weighted spherical harmonics expansion then becomes 

~ 2 ^ m = (\ - y) V^hn] - 7 E wl + [higher order ]. (25) 

where is defined as 

cfc?’ = I d£l- 2 Y? m (to) cos 6 - 2 Y Vm , (ft), (26) 

and has a non-zero values for £' = £ and (! — i =b 1 with m! = m given by 


2to r e+im 1 A - 1) A + 3) (l ~ m + 1 ) (l + m + 1 ) 
7^ + 1)’ ^ + l\/ (2^+1) (2^ + 3) 

(see also Appendix A of Ref. @3]). Because of the above result, we may consider C^" 1 ' = CpP n ,. 


(27) 


C. Improved extrapolation formula 


Comparing Eqs. 


and 


in the time domain, we have 


fj „ Al-W + 2)ft . lima *(*-l)(* + 2)(* + l) ~ 

-H-i.my' ^ JH- r ) . . llgmyt r ) -b o 2 H-imV' T ) 


2 r 


£ (^ + 1 ) r 


8r 2 


r ) T lrn Hi-lm(t r )^ 


= (77) [ r ^«m] - 7E C '^ / [ r ^4 R m'] + [higher order] 


la 

r 


( 28 ) 


£'m' 


Our improved extrapolation formula derived from the above equation is therefore 

- (l - 7) (-«.(*.<•> - (<, ~ 1 ([ < + 2) /(f,r)| 


+ 


(£- l)(£ + 2)(7 2 + £-4) 


8r 2 


J j dtdt[rip^(t,r )] 


2ia / (£ + 3)(^ — 1)(£ + m + 1)(£ — m + 1) 


(£+17 


(2f + l)(2f + 3) 


[rd t ipu+i m (t,r)\ - ^^- [np^ lm (t,r)] 
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2 i a j{£ + 2 ){i — 2)(U + m)(£ — m) 

W]J ( 2 £- 1 )( 2 £+ 1 ) 

+ [higher order]. 




(l ~ 2)(^ + 1) 

r 


WT-lm 



(29) 


The above formula ( [29] ) is our definitive equation for extrapolation of the waveform at finite radii to order 1/r 2 . It 
involves the first order correction in the mass (Schwarzschild-like) and spin (Kerr-like) of the sources, and corrects for 
the differences between the numerical and Kinnersley tetrads. 

On the other hand, we have proposed an extrapolation formula to order 1/r in [36] 


rV4 m 


r =oo 


(l _ ^ - (t 1] ^ + 2) J dlM®,(f,r)]) 

[ r 'p4p m '( t ’ r )] C em' ■ 

t'^t, m'=m 


(30) 


Since we did not take care of the difference between the spin-weighted spheroidal and spherical harmonics in the above 
equation, the spin correction is different between Eqs. (29) and (|30|). 


In the equations above r is the areal radius (Schwarzschild coordinate in the nonrotating case). In the standard 
numerical simulations we use Rnr that asymptotically behaves more like, R , the ‘isotropic’ radial coordinate, hence 
in Eq. (29) we typically use r = R( 1 + (M + a)/(2R)) (1 + (M — a)/( 2 R)). Alternatively, one could also compute 
directly the areal radius from the full numerical simulation via r = ^/A(R)/4tt, where A(R) is the measure surface 
area of the ‘sphere’ R = const. 


IV. ESTIMATION OF THE RADIATED ENERGY AND MOMENTA 


Using the improved extrapolation formula in Eq. (29), we derive extrapolation formulas for the radiated energy and 


momenta which are calculated from the Weyl scalar ?/q as 

iLJ dn J dtr ^ 

—E 

167r ^ 


dE 

dt 


dL z 

dt 


dPz_ 

dt 


1 


_sp 

167T 

J_c 

167r 


J dtripum\ r=l 
Jdtldc/) (/ dt rip 4 \ 
E m ( dt rip 4im \ r=c 


tm 

dCl COS 0 


/ dQ cos 9 [ dt r^ 41 

J 

Pm P'm' 


dt ril>Um\ r = c 


J J dtdtr^ | r=oc 


dtdt rijjum L_„ 


dt rip 


tm t'm‘ 


(31) 


Here, we focus only on the z component for the angular and linear momenta, and have used the normalization of ip 4 
as Eq. |l]). Cf™' is the same as in Eq. (26). 


Inserting Eq. (29) into the above expressions, we obtain the extrapolation formulas 

2 


dE 1 

rlt 1 (\ir E—' 


tm 
rd+lm 


-“)/ 


dt$e m 


+ 


~ 1)(^ + 2) 

2 r 2 


// 


dtdt^im 


-4a I V TO 5 ( $ 


£ +1 

(iLj; i r f, 

efa 167T \ 

+ (* -!)(* + 2) ( 


Ulm 

AM 
r 


/ 


—4a 


2r 2 

£ + 1 


dt&e 

j dt ^ j 

Jf***“) (Jff 


c: 


l—l m 

Am _c> i $ 


t—lm 


/ 


dtQgn 


dtdt&tr! 


dtdtdt&zm 


5? 




C: 


t—lm 
tm 


-5R 


(/A^'m) 

































dPz _ ]_ r l'm' 

dt 16tt ^ ^ lm 


irri I'm' 


(/ <lt ®' m ) (/■ 

/«><,„) (// dtdt&i'r 


Z) 


(£-£')(£ + £'+ 1 ) 

2 r 

2(^ - 1)^ + 1)(^ + 2) - ^ + !)£'(£' + 1) + 4 


4 r 2 


3? 


C; 


'l-lm 
lm o 


J J dtdt^em^j yj J dtdt$i'r 


■ 4a ( < /vr 9 (^+ im (/ 

2a / (2^ + 3) - (f - !){£' + 2))C*+ lm ^ 

- ITT - * 

(2(£ - 2)(£ + 1) - {£' - !)(£' + 2 ))Cf~ lm 


dt$ e+ dt$ e 

im 3 ? (( J dt$ t - im) (J dth' 


(32) 


where and we have ignored [higher order] terms described below Eq. (20). In order to simplify the 

expressions and to reduce the order of integration with respect to time, we have used the frequency domain analysis. 
In the expression of the radiated linear momentum, we take the sum over £' and m! as t = £, i =b 1 and m' = m. It 
should be noted that the (£, m) mode denotes the index of the spin-weighted spherical harmonics. There are order 
1/r corrections for the radiated energy and angular momentum that is different from Ref. [36] because of the tetrad 
difference. 


FULL NUMERICAL IMPLEMENTATION 


In order to evaluate the actual benefits of the analytic expression (30) for the extrapolation to infinity of gravitational 
waveforms extracted at a finite radii in a typical full numerical setting we consider the test case described in Table [T] 
We perform three sets of runs with increasing global resolution in the convergence regime and we extract waveforms 
at ten different radii, evenly separated as 1 /r. 


TABLE I. Initial data for our test case. The binary’s parameters were estimated using quasicircular orbits. 


Config. 

xi/M 

X2 /M 

P/M 

m\/M m v 2 /M 

S\/M 2 

s 2 /m 2 

mf/M 

mf/M 

Madm/M 

ai/rrii a^/m^ 

A_DU0.8 

-4.9832 

4.5267 

0.09905 

0.30178 0.30168 

-0.2 

0.2 

0.5 

0.5 

0.98951 

i 

o 

bo 

o 

oo 


In this work, we use a grid structure with 10 levels of refinement. The outer boundary was placed at 400M and 
for the medium resolution run the resolution was 4 M on the coarsest level and 1M in the wavezone. The finest level 
around each BH was as wide as twice the diameter of the relaxed horizon. We also performed a lower and higher 
resolution run with resolutions in the wavezone of M/0.88 and M/1.2. 

The simulation results will depend on the extraction radii as well as on the truncation errors due to finite resolution. 
Hence we consider different resolutions and extraction radii and extrapolations to null infinity. In this paper we 
used extraction radii up to R 0 \y S /M = 190 and locate the extraction radii equidistant in 1/R, with R 0 b s /M = 
75,80.4,86.7,94.0,102.6,113.0,125.7,141.7,162.3,190.0. 

We directly compared waveforms extracted with the characteristic method to our extrapolation formula, Eq. ©, 
in Ref. [16], Figs. 8-9, and to purely numerical extrapolations in Ref. (25). There we observed an excellent agreement 
with our analytic expression at first order in 1/r for the phase, as predicted by the error analysis of Eq. (11). The 
improvements in the amplitude are of higher order as shown in Eq. (10). In order to supplement those studies, here 


we focus on the integral expressions for the energy and momenta radiated at infinity. The results of such studies is 
displayed in Figs. [Tp] The radiated quantities are calculated using all modes up to i = 6 where the news and strain 
are calculated via the fixed frequency integration mm ■ 

In Fig. [l] we observe the computed radiated energy directly from the finite radii extraction that we denote as 
‘Raw’. The figure displays the different extraction radii, evenly distributed versus M/R 0 b s for the three finite- 
difference resolutions considered, denoted as Low, Medium and High. We provide a Richardson extrapolation to 
infinite resolution (3rd order) for each observer location value based on those three resolutions and also the value 
of the total radiated energy as inferred from the subtraction of the final horizon mass to the initial ADM mass of 
the system (denoted by the thick straight line). This measure of the final black hole mass, is very robust (at this 
scale) with increasing resolution and provides a very accurate measure as well as a consistency check of the extraction 
process. 
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Raw 


Perturbative 1/r 


Perturbative 1/r 2 





M/R n , 


M/Rn, 


M/Rn, 


FIG. 1. The energy radiated (adding up to i — 6) as a function of the observer location M/R 0 b s = 1/75, 1/80.4, 1/86.7, 1/94.0, 
1/102.6, 1/113.0, 1/125.7, 1/141.7, 1/162.3, 1/190.0 for the directly extracted waveform, labeled as ‘Raw’ (left) and for the 
analytically extrapolated waveform, labeled as ‘Perturbative 1/r’ (center) and ‘Perturbative 1/r 2 ’ (right). 


We observe that for the ‘Raw’ extraction increasing resolution (particularly for the closer to sources observers) 
brings the results further apart from the reference value inferred by the final horizon mass. To get consistency, one 
needs to first extrapolated to infinite observer location and then to infinite resolution. 

The second panel of Fig. [l] displays the same computation of the radiated energy, but after extrapolation of the 
waveforms via Eq. <§. We use the extrapolation at each observer location. We would expect that the dependence 
of the estimated energy radiated with the observer location is weaker since we are correcting for the 1/r behavior 
and only higher power dependencies should appear. We indeed observe flatter curves at all three finite-difference 
resolutions for this case compared to the ‘Raw’ extraction. The second feature is that at a single observer location the 
values converge towards the horizon value with increasing resolution. This is a desired feature, especially for a more 
demanding simulation where one only has access to accurate extraction in the intermediate zone between the sources 
and the radiation zone. The third panel shows the extrapolation carried to order 1/r 2 using Eq. (29) with a = 0 (in 
practice we did not see a strong dependence on a). Notably, in both cases, extrapolation to infinite resolution and 
infinite observer location leads to values within 0.1% of the correct value as inferred by the horizon measure. 

A similar behavior is observed in the computation of the angular momentum radiated as displayed in Fig. [2j For the 
first panel with the ‘Raw’ waveforms we see that increasing the finite-difference resolution leads to extrapolated values 
further apart from the horizon measure derived as the difference of the final spin of the black hole m to the initial total 
ADM angular momentum (denoted by the thick straight line). Using the perturbative 1/r and 1/r 2 extrapolations 
before the calculation of the angular momentum, as shown in the middle and left panels, leads to flatter curves with 
observer location and exhibits convergence toward the correct value with increasing finite-difference resolution. In 
both cases, the extrapolation to both infinite resolution and infinite observer location leads to predictions within 0.1% 
of the expected value. The importance of the extrapolation formula is that this can also be achieved with information 
from a single finite observer location. 
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FIG. 2. The Angular Momentum radiated (adding up to i — 6) as a function of the observer location M/R 0 b s = 1/75, 1/80.4, 
1/86.7, 1/94.0, 1/102.6, 1/113.0, 1/125.7, 1/141.7, 1/162.3, 1/190.0 for the directly extracted waveform, labeled as ‘Raw’ (left) 
and for the analytically extrapolated waveform, labeled as ‘Perturbative 1/r’ (center) and ‘Perturbative 1/r 2 ’ (right). 


Finally we also compute the linear momentum radiated by the system and display the results in Fig. [3] The first 
observation is that we do not have a very accurate measure on the final horizon for the recoil velocity to use as a 
reference value (although, see the work of Ref. @8]). However, based on the extrapolated values we estimate the recoil 
velocity to lie in the range 372 — 373km/s. We then observe that at a given finite value of the observer, particularly 
for those closer to the sources, the perturbative extrapolations values lie closer to the expected recoil. The curves also 
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look flatter indicating the internal consistency of the extrapolation process. 



FIG. 3. The Linear Momentum radiated (adding up to i — 6) as a function of the observer location M/R 0 \> s = 1/75, 1/80.4, 
1/86.7, 1/94.0, 1/102.6, 1/113.0, 1/125.7, 1/141.7, 1/162.3, 1/190.0 for the directly extracted waveform, labeled as ‘Raw’ (left) 
and for the analytically extrapolated waveform, labeled as ‘Perturbative 1/r’ (center) and ‘Perturbative 1/r 2 ’ (right). 


In order to produce a reference waveform that we may consider the best extrapolation and hence approximation 
to the exact waveform in Fig. [4j we took the highest resolution run and used the ten extraction radii we have to 
extrapolate the waveform in time using a 2nd order fitting polynomial in 1/r. We extrapolated the amplitude and 
phase after shifting the time by the tortoise radius for each extraction radius. We then can compare the amplitude and 
phase of this extrapolated waveform to a finite radius waveform (i7 0 bs = 190M, our largest extraction radius), and to 
the waveforms produced by using the 1/r and 1/r 2 order perturbative extrapolations (without the terms depending 
on the spin). The results are displayed in Fig. p] which shows the benefits of using our formulas to approximate the 
waveform phase and amplitude at infinity. Note that given the different dependence of the phase correction (1/r as 
shown in Eq. @» and the amplitude correction (1/r 2 as shown in Eq. ([To|) the phase and amplitude show further 
improvements by including the second order corrections. This is more explicitly displayed in Fig. [6j that summarized 
the averaged differences. 



t/M 



t/M 



t/M 



t/M 




t/M t/M 


FIG. 4. ‘cx)’ is the extrapolated NR waveform (phase and amplitude) at highest resolution using 10 radii equally spaced between 
75 and 190M. We compare the above best extrapolated waveform, ‘oo’, with that directly extracted at R 0 bs = 190M, 1st 
order, (1/r), perturbative extrapolation (PE), and 2nd order,(1/r 2 ), PE. In all cases the mode (£, m) = (2,2) is displayed. The 
first row is the Weyl Scalar rMT 22 and the second row is the gravitational strain rh 22 /M. 
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t/M t/M 

FIG. 5. Comparison of the best extrapolated waveform, ‘oo’, with that directly extracted at R 0 b s = 190M, 1st order, (1 /r), 
perturbative extrapolation (PE), and 2nd order, (1/r 2 ), PE, calculating the difference abs(‘oo’-waveform) for the phase and 
amplitude. The mode (£, mn) = (2, 2) is displayed. 




^obs^ ^obs^ 


FIG. 6. (£,m) = (2,2) mode of the raw (NR) and first (EX1) and second (EX2) perturbative extrapolated waveforms as a 
function of radius versus the waveform extrapolated to infinity. Displayed is the mean of the % difference between perturbative 
and extrapolated for each radii for the amplitude and phase for the times between 400 M and 800M. 


VI. DISCUSSION 


In this paper we describe a procedure for extrapolating the waveform at finite radius to infinity as a power series 
in 1/r. We provided the complete 1/r correction to the waveform in Eq. ( |30| ), including the spin terms in the 
background. We have also found it important to include the leading terms in 1/r 2 in the extrapolation formula 
as given in Eq. ( [29] ) . We have tested the extrapolation formula’s properties in a typical full numerical simulation 
of a black-hole binary, where we can verify the behavior of the extrapolation with different observer locations and 
finite-difference resolutions. In numerical simulations where the typical wavelength is relatively small compared to 
the boundaries of the simulation, the perturbative extraction provides at least a way of verifying the accuracy and 
consistency of the waveforms and radiative quantities such as the total energy and linear and angular momenta. In 
situations where it is extremely costly or inaccurate to extract at distances of two gravitational wavelengths from 
the sources (rule of thumb for the radiation zone), this method provides a crucial technique to evaluate waveforms 
and radiated quantities. In particular, we have seen that it is only the extrapolated waveform that converges with 
increasing resolution to the correct values and that extrapolation to infinite resolution of a finite extraction waveform 
can lead to a worse approximation. Although, for far enough location observer and resolution these two extrapolation 
processes eventually tend to commute. 


The second order correction provided in Eq. (29) could be useful in situations where we have extended sources 
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or one needs extreme resolutions near the sources and the simulation grid cannot reach the radiation zone. It also 
provides an independent way to estimate the errors of the first order extrapolation formula Eq. (30) by looking at the 
differences produced by these two extrapolation formulas. 
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Appendix A: Second order correction with ^o m 


In Sub-Sec. 

the identities in the Teukolsky formalism, and use the notation as given in Ref. [49] . In the Schwarzschild background 
with mass M, the Weyl scalar ^4 and V’o have following relations. 

r 4 ^4 = ^2 r4 ^ 2 (j dt ~ 5r ) r4 ^ 2 ^ ’ V’o = g (Yf + 12 Mdt^j , (Al) 

where / = 1 — 2 M/r and 9 = —(do — scot# + icscOd$) for the spin-s weighted spherical harmonics [50] . 4> denotes 
a Hertz potential. 

Here, since we are interested in the leading asymptotic behavior for large r, the equation for ^4 is approximated as 

r ^ 4 = , (A2) 

where T = (d t — 9*). In the above equation, the left hand side is written with respect to the retarded time t — r*. 
Therefore, using T -1 = ( 1 / 2 ) / dt, we have 

dt dt dt dt (r?/q) = . (A3) 

For -0o, we ignore the term proportional to M in order not to introduce the complex conjugation of 4^ m , and focus 

4 _ _ 

on the term 9 4/ is a spin-(— 2 ) function, i.e., 

Y[ f f [ dtdtdtdt(r^i m )-2Yi m = ^r 5 Y^ £m -2 Y im- (A4) 

im ** J J im 



IIB, the formula has a term, ^4 integrated twice in time. In order to remove this term, we may use 


The operator 9 on 4t gives 


- l W + W + 2)^ m 2 Y lm , 

im 


(A5) 


where there is no change in the (t, r) dependence because 9 acts only on the angular variables. Although there may 
be a relation between 12 M9 t ^ and the term proportional to M in Eq. (A 8 ) below because both of the numerical 
factors are 3/2, we simply ignore it here in order not to introduce the complex conjugation of 4^ m . This means that 
we consider an approximation, 


Y V>o m 2 Y lm = \Y^~ +W + 2)^ m 2 Y tm . 

im im 


Combining Eqs. (A4) and (A 6 ), we have for each (£, m ) mode 


1 

8 


(t- 


l)^ + l)(* + 2) //// dt dt dt dt (ri/jf n ) 


1 5 / im 
2 V % ' 


(A 6 ) 


(A7) 
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in the large r limit and the above approximation. Therefore, Eq. (16) is rewritten as 


i.tm I 


V Ip A I 


rV 4 m (t,r) - J dt[r^ 4 m (t,r)} 


+ 


(^ 2 + £ - 4) 
2 £(£ + 1 ) r 2 


3 M 

#[r 5 ^ m (*,r)] -j^Jdt [r^ m (t,r)] + 0(l/r 3 ). (A 8 ) 


Here, we have used V ; o m (^ r ) extracted at a finite radius because the error due to the use of fi nite extraction radii 
becomes higher order in the large r expansion. Since we have used an approximation to derive Eq. (A 6 ), for consistency, 
the M-dependent term should not be kept any more, i.e., 


(£-!)(£+ 2 ) 


rVCUoo =r^(t,r)~ 2r 

(£ 2 + £ — 4) 2 r 5 . Pm 


J dt[riPi m (t,r )] 


+ 


2 £[£ + 1 ) r 2 


dtv ro m (t,r)\ + [higher order] . 


(A9) 


This derivation is in an ideal situation where we have assumed that there is no contribution from the other Weyl 
scalars, the peeling theorem applies, and we have used a low frequency Muj approximation. 
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